{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Estimating current cases in Wuhan\n",
    "\n",
    "#### Author: Yiran Jing\n",
    "#### Date: Jan 2020\n",
    "\n",
    "Based on cases detected outside mainland China. **We use Observation on Jan 22**. Since on 2 am Jan 23, Wuhan shut down.\n",
    "\n",
    "\n",
    "## Main Conclusion:\n",
    "- There are at least **4600** 95% CI(2100, 8550) cases in Wuhan, until Jan 19. (There is time deley between suspected and confirmed)\n",
    "- Based on 29 confirmed overseas cases in Jan 26, There are at least more than **16600** cases 95% CI (11310, 23480) before Jan 23.\n",
    "- Based on 67 confirmed overseas cases in Jan 29, There are at least more than **38500** cases 95% CI (30000, 48470) before Jan 23.\n",
    "- **Commuting flows** has significant impact on `2019-nCov epidemic growth rate`\n",
    "***\n",
    "\n",
    "## Statistical Modelling:\n",
    "\n",
    "### Sensitivity analysis: \n",
    "Sensitivity analysis to estimate current cases in wuhan based on 3 scenario\n",
    "1. Baseline\n",
    "     - **8** overseas confirmed cases until 22 Jan.\n",
    "     - 10 mean time to detection\n",
    "     - 19 million airportCatchment\n",
    "2. Smaller catchment:\n",
    "     - airportCatchment = wuhan_population = 11 million\n",
    "3. Shorter detection window:\n",
    "     - 8 mean time to detection\n",
    "\n",
    "### Profile likelihood CI\n",
    "In general, the confidence interval based on the standard error strongly depends on the assumption of normality for the estimator, which is something we cannot guarantee in this case. The \"profile likelihood confidence interval\" provides an alternative.\n",
    "\n",
    "We defined a Binomial likelihood for the number of exported cases and used this function to find the MLE of the number of cases in Wuhan, using the profile likelihood approach to identify the 95% CI around the MLE. The lower and upper bounds of the 95%CI are those values by which the log-likelihood difference from the maximum log-likelihood is 1.92 (95%-percentile of the Chi-square(1) distribution).\n",
    "\n",
    "\n",
    "***\n",
    "\n",
    "\n",
    "### Reference: \n",
    "1. [Estimating the potential total number of novel Coronavirus cases in Wuhan City, China (Jan 21 2020)](https://www.imperial.ac.uk/media/imperial-college/medicine/sph/ide/gida-fellowships/2019-nCoV-outbreak-report-22-01-2020.pdf)\n",
    "2. [Confidence intervals by the profile likelihood method](http://people.upei.ca/hstryhn/stryhn208.pdf)\n",
    "3. [Binomial profile likelihood CI ](https://personal.psu.edu/abs12/stat504/Lecture/lec3_4up.pdf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from dataclasses import dataclass\n",
    "from typing import List\n",
    "from scipy.stats import nbinom, t\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import stats\n",
    "import random\n",
    "import math\n",
    "from math import lgamma\n",
    "from scipy.optimize import minimize\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "from help_function_model1 import * "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2019-01-21 \n",
      "Scenario: Baseline\n",
      "Exported number of confirmed cases: 7\n",
      "Detection Window: 10\n",
      "Time from onset of symptoms to detection: 4\n",
      "Daily International passengers travelling out of Wuhan: 3301\n",
      "Effective catchment population of Wuhan Interntional Airport: 19000000\n",
      "\n",
      "Estimated number of cases in Wuhan: 4029.0820963344436\n",
      "Estimated 95% confidence interval: (1730, 7780)\n",
      "\n",
      "------------------------------------------------------------------------------------\n",
      "\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAEICAYAAACj2qi6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3df5xcdX3v8deb/AIJJpCkEAmXhEvEBqtItwitVSsIgbbG2tx7E70KCqUq1LbYH3DtVUuxj3Jvb/Fq+VEK3CIFAqUBghdFr9ASRIHF8MOAS9aEmGAWlvyChPza7Of+cb4bJsPM7tndmZ05M+/n4zGPPfM953zP9zs7M5/z/XHOKCIwMzPL44BGF8DMzIrDQcPMzHJz0DAzs9wcNMzMLDcHDTMzy81Bw8zMcnPQaFKSDpJ0j6Stkv5lBPt/RtKLkrZJmlaPMjaSpJWS3t/ocjSKpGsk/fcmKMf7Ja1vdDls7DhoNK+FwOHAtIj4T8PZUdIE4O+A0yNickRsrEcBGykijo+If2t0ORolIj4dEX813P0kLZb0bFnad6ukXTzactaSpD+V9GNJr0paI+lPy9bPlvSApNck/UTSaSXr3i7pPkkvS3rDxWmSflHS/ekkrVvS7wxRlsMk3Slpu6S1kj5asm6mpGWSfi4pJM0eIq/flPSQpC2SeiRdJ+mQkvX/WdLDqV7/NuQLVWcOGs3raOC5iOgbwb6HAwcCK2tbpP1JGl/P/EeqWcvVJB4E3iZpBux7rd4JHFSWdkratpkI+ARwKDAfuFDSopL1twIrgGnAF4A7BuoE7AFuB859Q6ZZfe8GvgkcBpwP/LOktw5SliuB3WSftY8BV0s6Pq3rB74N/G7Oek0BLgPeAvwicCTwP0vWbwK+CvxNzvzqKyL8qNEDCOBzwGrgZbJ//AFp3TnA94ErgC1pm19N6euAl4Cz07Z/SfaG3ANsA86tcKxJZG+kn6fHV1PaW4HtqSzbgPurlPU9wMOpLOuAc1L6FOAbQC+wFviLKnXYSPZGPyBtszbV4RvAlLT97FSOs4GfpdfkCyVlOAn4QSrDBuDvgYlp3dXA35aV+W7gorT8PHBaWv4ycAfwz8ArwHnAPwGXlez7fmB9yfM/B14AXgW6gFOrvE4HAf8r1W8r8BBwUFr3L0BPSn8QOL5kv7OAZ1L+LwB/UrLut4AnUr0fBt4xgnLtq99A3YDPp//BBuCTg7xPfwr8bsn/4AHgxrK014AJJe/rY0dybOA3yb7IXyF7n325ZN2g748cn7evAV9Py28FdgGHlKxfDny6bJ9jgShLezvZZ0Ulad8B/qrKcQ8m+3y+tSTtJuBvyrYbn+o3e5jfIx8Bnq6Qfh7wb7X+3hruwy2N2vsdoAM4EVgAfKpk3buBp8jOhG4BlgC/QvZG/q/A30uaHBFfAv4auC2y7qXrKxznC8DJwAlkZ4onAX8REc8BA2c8UyPiA+U7Sjoa+BbwdWBGyuOJtPrrZIHjGOB9ZGd2nyyrw2qyM6yvkAWSc4DfSPtMJvvyL/Ue4DjgVOCLkn4xpe8F/hiYTnZmeyrw2bTuVuC/SFIq86HA6ek1q2QBWeCYCtxcZZuB+h8HXAj8SkQcApxBFoQq+Vvgl8kC/GHAn5GdSUL2Gs4FfgH4Udlxrwd+P+X/duD+dOx3ATcAv0/2PvgHYJmkScMsV7kjyP5vR5KdTV+ZXrNKHgTem5bfS/bl+lBZ2g8jYk8Njr2d7D00lSyAfEbSh8v2r/b+qCq9L36d11vTxwOrI+LVks2e5PXPwnCJ7P9WyVuBvvRZq8Wxyr2XOvcSjEqjo1YrPcjOKuaXPP8s8L20fA6wqmTdL6XtDy9J2wickJa/DPzzIMf6KXBWyfMzgOfT8uyU9/gq+14C3FkhfRzZGdS8krTfJ53dpDr8rGyf7wGfLXl+HFkLaXxJOWaVrH8UWFSlXH80UC6yD+3PgPem579HSauJN7Y0HizL65+o0tIgC9IvAaeRzqarlOcAYAfwzhz/+6mprgOtrJ+l1+7NZdtdTdkZLFmL4n15y1Vev1S3HaX/75TPyVX2PQdYkZbvBj4IvK0s7Utl7+vBWhrDOfZXgSvK3qe53h9l+fwl2Rf1pPT842SBrnSbrwD/VJZWqaUxgexE6M/S8ulkn4P7qhz714GesrTfo6wVwAhaGul/sZmSVkzJOrc0WtS6kuW1ZP2UA14sWd4BEBHlaZNzHuctKf9qxxrMUWRBp9x0sg9Neb5Hljxfx/4qlWM8WUtkQE/J8mukOkp6q6RvpsG/V8haV9MhfaqzVsXitN9HGbwFUV6uqiKimyxAfRl4SdISSZVeu+lkY0NveK0kjZP0N5J+msr+fMk+kPVnnwWslfTvkk5J6UcDn0+DnlskbSH7f7xlGOWqZGPsP/6173Wu4EHgHak1cDLwg4j4CTAzpb2H4Y1nVD22pHenweleSVuBT/P6azSg4vujGkkXkrVefjMidqXkbcCbyzZ9M1k336Aia1F9mKwl1EPW1XY7Wbcbkr6VZiFuk/Sx0RxL0q+X5LWybN3JZD0QC2P/VkxTcdCovaNKlv8D2XhDPfyc7AtoJMdaB/zHCukvk7USyvN9oeR5+cyTSuXoY/8AWc3VwE+AuRHxZuC/kbUwBtwKLEzdae8G/nWQvMrLtR14U8nzI/bbOOKWiHhPKnsAl1fI82VgJ5Vfq4+SdYmdRtY1MzulK+X/WEQsIOu6uovsSwiy1/4rETG15PGmiLh1GOUalYhYTfZ/O5+s5bgtrfpBSpsM/LBkl9cY5LUcwi3AMuCoiJgCXMP+/+NhkfQp4GKysZ7Sqb4rgWNKZx2Rddvm6uaJiKci4n0RMS0iziDran00rTszsm7iyRFxM/AcMF7S3OEeKyKWl+S1rzsrdVsuAz4VEd/LU+ZGcdCovT+VdKiko4A/BG6r03FuBf5C0gxJ04Evkg0E53EzcFqayjde0jRJJ0TEXrIvt69IOiR9WV80RL63An8saY6kybw+FpNn1tchZAOk2yS9DfhM6cqIWEH2xX0dWVfBlpz1g2yM5qw0NfIIsjN4IBvTkPQBSZPIgsIOXh+nKD1+P9n4w99JektqXZyS9juEbOB1I9kX6l+X5D9R0sckTUlnsa+U5P+PwKfTGbgkHZymXB6St1w1spzsf7u8JO2hlNYZETtK0p8APprqP5+sKy2vQ4BNEbFT0klkwXZE0ln+XwMfTIFvn3Rm/gTwJUkHKpsy+w7SiUZ6rQ8EJqbnB6bXeSDvd6S0N0n6E2AmWTfcG0TEdmApcGn6//0a2QnETSX5HUg2MQVgUnperV5vJ5tt9QcRcU+F9ePS/uOBA1I5J1R/perLQaP27gYeJ3sD/1+yAdF6uAzoJBtYf5psIPayPDtGxM/Iuk4+Tzad7wmyMyWAPyA7S19N9iVyC9kXZzU3kH1YHgTWkH3Z/UHOOvwJ2ZfIq2RfppUC7C1kZ/O35MxzwE1kfd7Pk82EKc17Etn0xZfJuiN+gWycp1oZnwYeI3utLif73HyDrCvuBbJZUj8s2+/jwPOp6+rTZNMyiYhOsv7vvyfru+4mG2MYbrlG699T/g+VpC1PaeVdU38I/DbZbK+PkbWc8vos2Zfrq2QnNrcPsf1gLiObPPBYSRfPNSXrF5FNQtlM9joujIjetO5osiA80BrYQTaWNODjZDO/XiIbkP9gSddXtXodlLa/FfhMRJS2NHaQdWNB1preQXWfJ5uQcn2VrquPp/2vJhtP2UH2eWkIpQEWqwFlFw3NTX3TZmYtxy0NMzPLzUHDzMxyc/eUmZnl5paGmZnl1tI3dps+fXrMnj270cUwMyuUxx9//OWImFFpXUsHjdmzZ9PZ2dnoYpiZFYqktdXWuXvKzMxyc9AwM7PcHDTMzCw3Bw0zM8vNQcPMzHJz0DAzs9wcNMzMLDcHDTMzy62lL+6z5rZq0zae7h3yFzL3mTzhAE4/5vChNzSzunHQsDFzV9eGUf0M3bY9/Szt2rBf2keOmzm6QpnZsDhoWF09sLaXzTvz/PLryJQGEQcQs/pz0LC6uLe7h517x/a2+wMB5IiDJ/Krs6aN6bHN2oWDhtXUip4trNk62M8h11/P9t0s7drg4GFWB549ZTWztGtDwwNGqYHgYWa145aGjVq9xy1GayBweMzDbPTc0rBRWdq1oakDRqmlXRt4eP3GRhfDrNAcNGzEitj14y4rs9Fx0LBh631tV+G/eJd2baD3tV2NLoZZ4XhMw4alGWZH1crydZuYcAD89lyPdZjl5ZaG5fbw+o0tEzAG7OkvZjebWaM4aFguD6/fSM/23Y0uRt04cJjl46BhQ2r1gDHA4xxmQ3PQsEGt2rStLQLGgOXrNvHA2t5GF8OsaTloWFW9r+0a1q3LW8XmnX3c293T6GKYNSUHDatq+bpNjS5Cw+zcG9yzyuMcZuUcNKwiDwxnM6vu8utgtp9cQUPSfEldkrolXVxh/SRJt6X1j0iaXbLukpTeJemMofKUdHNK/7GkGyRNSOnvl7RV0hPp8cXRVNyqu9NflPv048BhVmrIoCFpHHAlcCYwD1gsaV7ZZucCmyPiWOAK4PK07zxgEXA8MB+4StK4IfK8GXgb8EvAQcB5JcdZHhEnpMelI6mwDe6eVRsY21/BaH79OJCaDcjT0jgJ6I6I1RGxG1gCLCjbZgFwY1q+AzhVklL6kojYFRFrgO6UX9U8I+LeSIBHgVmjq6LltaJnC3tG83usLSxwl50Z5AsaRwLrSp6vT2kVt4mIPmArMG2QfYfMM3VLfRz4dknyKZKelPQtScdXKqyk8yV1Surs7fXUyeFotau968EtDmt3zTwQfhXwYEQsT89/BBwdEe8Evg7cVWmniLg2IjoiomPGjBljVNTi81l0PoEDh7W3PEHjBeCokuezUlrFbSSNB6YAGwfZd9A8JX0JmAFcNJAWEa9ExLa0fC8wQdL0HOW3IXhq6fAEHhy39pUnaDwGzJU0R9JEsoHtZWXbLAPOTssLgfvTmMQyYFGaXTUHmEs2TlE1T0nnAWcAiyNiXw+7pCPSOAmSTkpl9y/qjNKqTds8jjEC/TjYWnsa8tboEdEn6ULgPmAccENErJR0KdAZEcuA64GbJHUDm8iCAGm724FngD7ggojYC1Apz3TIa4C1wA9SjFiaZkotBD4jqQ/YASxKgclGoR2v+K6VPf1wb3cPZx17RKOLYjZm1Mrfux0dHdHZ2dnoYjStu7o24EbG6E2ecACnH3N4o4thVjOSHo+Ijkrrmnkg3Oro4fUbHTBqZNuefv/2uLUNB4021U53rh0LPdt3s2rTtkYXw6zuHDTakGf+1IfHh6wdOGi0mRU9W9wtVUe+3sVanYNGm/FV3/XnwGGtzEGjjfi6grHjwGGtykGjjfgivrF193MOHNZ6HDTahM98x97egO+sfrHRxTCrKQeNNvDAWt/tt1G27elnRc+WRhfDrGYcNNrA5p19jS5CW/PkA2slDhotzv3qzcHdg9YqHDRa3N7WvbVY4ThwWCtw0Ghh/rGg5uOr8a3oHDRa1IqeLbiR0Xz68cQEKzYHjRblwdfmtXlnn29uaIXloNGCfCbb/HxzQysqB40W5Cm2xeCBcSsiB40Wc293T6OLYMPggXErGgeNFrPTc2wLpR/fasSKxUGjhfhCvmLyrUasSBw0WogbGcXl2W5WFA4aLcKtjOLzwLgVgYNGi3ArozV4YNyanYNGC/AXTevwwLg1OweNgut9bRf+Qb7Wsm1PP72v7Wp0McwqctAouOXrNjW6CFYH/r9as8oVNCTNl9QlqVvSxRXWT5J0W1r/iKTZJesuSeldks4YKk9JN6f0H0u6QdKElC5JX0vbPyXpxNFUvBX4bLS1eWDcmtGQQUPSOOBK4ExgHrBY0ryyzc4FNkfEscAVwOVp33nAIuB4YD5wlaRxQ+R5M/A24JeAg4DzUvqZwNz0OB+4eiQVbiUP+Wy05fkKf2s2eVoaJwHdEbE6InYDS4AFZdssAG5My3cAp0pSSl8SEbsiYg3QnfKrmmdE3BsJ8Cgwq+QY30irfghMlTRzhPUuvN7XdvnW521g597whX/WVPIEjSOBdSXP16e0ittERB+wFZg2yL5D5pm6pT4OfHsY5UDS+ZI6JXX29rbu3V7d590+fOGfNZNmHgi/CngwIpYPZ6eIuDYiOiKiY8aMGXUqWmN5LKP9eHzDmkWeoPECcFTJ81kpreI2ksYDU4CNg+w7aJ6SvgTMAC4aZjnaglsZ7clX/VszyBM0HgPmSpojaSLZwPaysm2WAWen5YXA/WlMYhmwKM2umkM2iP3oYHlKOg84A1gcEf1lx/hEmkV1MrA1IvwpsraxN/wDW9Z444faICL6JF0I3AeMA26IiJWSLgU6I2IZcD1wk6RuYBNZECBtdzvwDNAHXBARewEq5ZkOeQ2wFvhBNpbO0oi4FLgXOItsMP014JO1eAGKxld/t7fNO/vofW0XM940qdFFsTalrEHQmjo6OqKzs7PRxagp920bwEeOa9uJgzYGJD0eER2V1jXzQLiVcZ+2DbjTJw/WIA4aBeI72dqAwBf+WWM4aBSEWxlWzhf+WSM4aBSEWxlWiS/8s7HmoFEA7oawwXhyhI0lB40C2Olmhg3hnlUOHDY2HDSanH/FzfLY0w8Pr9/Y6GJYG3DQaHLb9vh3+Syfnu27G10EawMOGk3Mt4yw4fL4htWbg0YT27yzr9FFsALyrWasnhw0mpTn39tI9eNWqtWPg0aT8vx7G42BGxua1ZqDRhPyh91qwb+7YvXgoNGEHvKH3WrEA+NWaw4aTciX8lkt+Y4CVksOGk3GV/ZarfnGhlZLDhpNxtfyWT14YoXVioNGE/EtQ6yePL5hteCg0UR8yxCrN1/4Z6PloNEkfDGWjYV+3KK10XHQaBK+ZYiNlW17+lm1aVuji2EF5aDRBPwBtrH2dO+rjS6CFZSDRhPwB9gawQPjNhIOGmZt7O7nHDhseBw0GswfWmukveFJGDY8DhoN5p//tkbzHXFtOBw0Gsj3BLJm4TviWl65goak+ZK6JHVLurjC+kmSbkvrH5E0u2TdJSm9S9IZQ+Up6cKUFpKml6S/X9JWSU+kxxdHWulmsdPNDGsiHhi3PIYMGpLGAVcCZwLzgMWS5pVtdi6wOSKOBa4ALk/7zgMWAccD84GrJI0bIs/vA6cBaysUZ3lEnJAelw6vqs3l4fUbG10EszfwGJsNJU9L4ySgOyJWR8RuYAmwoGybBcCNafkO4FRJSulLImJXRKwBulN+VfOMiBUR8fwo69X0erbvbnQRzN7AA+M2lDxB40hgXcnz9Smt4jYR0QdsBaYNsm+ePCs5RdKTkr4l6fhKG0g6X1KnpM7e3uZ88/tiPmtmm3f2+T1qVRVpIPxHwNER8U7g68BdlTaKiGsjoiMiOmbMmDGmBczLF/NZs/N71KrJEzReAI4qeT4rpVXcRtJ4YAqwcZB98+S5n4h4JSK2peV7gQmlA+VmVlseGLdK8gSNx4C5kuZImkg2sL2sbJtlwNlpeSFwf0RESl+UZlfNAeYCj+bMcz+SjkjjJEg6KZW9cKPJ/mU+KxLfSt3KDRk00hjFhcB9wLPA7RGxUtKlkj6UNrsemCapG7gIuDjtuxK4HXgG+DZwQUTsrZYngKTPSVpP1vp4StJ16RgLgR9LehL4GrAoBaZC8U9mWJH04+uJbH8q4Pdubh0dHdHZ2dnoYuzzwNpe3wLdCmnOlIN41xFTG10MGyOSHo+IjkrrijQQXngOGFZU/o1xG+CgMUZW9GxpdBHMRsUD4wYOGmPGZ2rWChw4zEFjDPgOotZKPKOqvTlojIGHfAdRayH9eOp4O3PQGAOtOz/N2tWeft+jql05aNSZ57hbq/I9qtqTg0ad+TczrJX5HlXtx0Gjjtx8t3bgGVXtxUGjjnwxn7ULB4724aBRJ+7rtXZzpwNHW3DQqBP39Vq7Cfxzse3AQcPMamZvwHdWv9joYlgdOWjUgS98sna2bU+/J4G0MAeNOvBvZli727yzzzfpbFEOGjXmprlZZs3WHb7vWgty0KixbW5mmO2z3PddazkOGjX08PrC/WS5Wd35Go7W4qBRQz3bdze6CGZNyYGjdTho1Igv5jMbnANHa3DQqBFfzGc2NAeO4nPQMLMx5cBRbA4aNeBbJ5gNj+9TVVwOGjXgn8wwG57AgaOoHDRGyRfzmY2MA0cxOWiMki/mMxu5AO5y4CiUXEFD0nxJXZK6JV1cYf0kSbel9Y9Iml2y7pKU3iXpjKHylHRhSgtJ00vSJelrad1Tkk4caaVrxTdlMxu9fhw4imTIoCFpHHAlcCYwD1gsaV7ZZucCmyPiWOAK4PK07zxgEXA8MB+4StK4IfL8PnAasLbsGGcCc9PjfODq4VW19vzLfGa10Y8nlBRFnpbGSUB3RKyOiN3AEmBB2TYLgBvT8h3AqZKU0pdExK6IWAN0p/yq5hkRKyLi+QrlWAB8IzI/BKZKmjmcytaSbxliVlt7w4GjCPIEjSOBdSXP16e0ittERB+wFZg2yL558hxJOZB0vqROSZ29vfXrPvItQ8xqb2+4q6rZtdxAeERcGxEdEdExY8aMuhzDtwwxqx+PcTS3PEHjBeCokuezUlrFbSSNB6YAGwfZN0+eIynHmPAtQ8zqqx9Px21WeYLGY8BcSXMkTSQb2F5Wts0y4Oy0vBC4PyIipS9Ks6vmkA1iP5ozz3LLgE+kWVQnA1sjYszfVf5RGbOxEfiWI81oyKCRxiguBO4DngVuj4iVki6V9KG02fXANEndwEXAxWnflcDtwDPAt4ELImJvtTwBJH1O0nqylsRTkq5Lx7gXWE02mP6PwGdHXfsReMg/KmM2phw4mouyBkFr6ujoiM7Ozprm6TewWWN85LiGTZZsO5Iej4iOSutabiC8njwd0KxxfMLWHBw0hsE3JjRrLAeOxnPQyOne7p5GF8HMcOBoNAeNnHa6mWHWNJZ2bfD1Ug3ioJGDWxlmzefp3ld909AGcNDIwa0Ms+a0eWcf96xyd9VYctAYglsZZs1tT79vOzKWHDSG4FaGWfPrxwPkY8VBYxD+KVezYnHgqD8HjUH4p1zNisczq+rLQaMKtzLMisszq+rHQaMKtzLMim3zzj4PkNeBg0YFbmWYtQYPkNeeg0YFbmWYtZalXRv8Wzg14qBRxgNoZq1p+bpNvu6qBhw0yvx8285GF8HM6mTn3nB31Sg5aJhZ23HgGDkHDTNrS0u7NnjSywg4aJhZ29q2p9+tjmFy0DCztufZVfk5aJiZkc2u8sWAQ3PQMDNLfDHg0Bw0zMzKLO3a4B93qsJBw8ysgj39bnVU4qBhZjaIpV0buPs5B48BDhpmZkPYG251DMgVNCTNl9QlqVvSxRXWT5J0W1r/iKTZJesuSeldks4YKk9Jc1Ie3SnPiSn9HEm9kp5Ij/NGU3Ezs+Fa2rWBO9s8eAwZNCSNA64EzgTmAYslzSvb7Fxgc0QcC1wBXJ72nQcsAo4H5gNXSRo3RJ6XA1ekvDanvAfcFhEnpMd1I6qxmdkoBFnweHj9xkYXpSHytDROArojYnVE7AaWAAvKtlkA3JiW7wBOlaSUviQidkXEGqA75Vcxz7TPB1IepDw/PPLqmZnVR8/23W3ZZZUnaBwJrCt5vj6lVdwmIvqArcC0Qfatlj4N2JLyqHSs35X0lKQ7JB2Vo+xmZnXVbl1WRRoIvweYHRHvAL7L6y2b/Ug6X1KnpM7eXv9GsJnV30CXVTv8LnmeoPECUHpWPyulVdxG0nhgCrBxkH2rpW8EpqY89jtWRGyMiIGbw1wH/HKlwkbEtRHREREdM2bMyFE9M7Pa2Lyzr+XvY5UnaDwGzE2zmiaSDWwvK9tmGXB2Wl4I3B8RkdIXpdlVc4C5wKPV8kz7PJDyIOV5N4CkmSXH+xDw7PCqamY2Npav29Sy4x3jh9ogIvokXQjcB4wDboiIlZIuBTojYhlwPXCTpG5gE1kQIG13O/AM0AdcEBF7ASrlmQ7558ASSZcBK1LeAJ+T9KGUzybgnFHX3sysjgYCx0eOmznElsWh7OS+NXV0dERnZ+ew9vn3n73Mxh176lQiM2tnRQkekh6PiI5K64o0EG5mVmhLuzYUvtvKQcPMbIwVOXg4aJiZNUgRg4eDhplZgxUpeAw5e8rMzMbGQOCYM+Ug3nXE1AaXpjIHDTOzJrNm6w7WbN3BhAPgt+c214wrd0+ZmTWpgV8PbKauK7c0zMwKYCBwHHrgeH7j6MbdIsktDTOzAhm4v1WjWh9uaZiZFdRA4DgA+PAYXW3uoGFmVnD9vB5AJk84gNOPObxux3LQMDNrIdv29O8LIEccPJFfnTWtpvl7TMPMrEX1bN/Nd1a/WNM8HTTMzFrYtj39Nc3PQcPMzHJz0DAzs9wcNMzMLDcHDTMzy81Bw8zMcnPQMDOz3Bw0zMwsNwcNMzPLzUHDzMxyc9AwM7PcHDTMzCw3Bw0zM8vNQcPMzHLLFTQkzZfUJalb0sUV1k+SdFta/4ik2SXrLknpXZLOGCpPSXNSHt0pz4lDHcPMzMbGkEFD0jjgSuBMYB6wWNK8ss3OBTZHxLHAFcDlad95wCLgeGA+cJWkcUPkeTlwRcprc8q76jHMzGzs5GlpnAR0R8TqiNgNLAEWlG2zALgxLd8BnCpJKX1JROyKiDVAd8qvYp5pnw+kPEh5fniIY5iZ2RjJEzSOBNaVPF+f0ipuExF9wFZg2iD7VkufBmxJeZQfq9ox9iPpfEmdkjp7e3tzVG9/b5l84LD3MTNrVpMn1HbouuV+IzwirgWuBejo6Ijh7j/3sMnMPWxyzctlZtYK8oSgF4CjSp7PSmkVt5E0HpgCbBxk32rpG4GpKY/yY1U7hpmZjZE8QeMxYG6a1TSRbGB7Wdk2y4Cz0/JC4P6IiJS+KM18mgPMBR6tlmfa54GUBynPu4c4hpmZjZEhu6ciok/ShcB9wDjghohYKelSoDMilgHXAzdJ6gY2kQUB0na3A88AfcAFEbEXoFKe6ZB/DiyRdBmwIuVNtWOYmdnYUSufrHd0dERnZ2eji2FmViiSHo+IjkrrfEW4mZnl5qBhZma5OWiYmdcGef8AAASVSURBVFluDhpmZpZbSw+ES+oF1ja6HMl04OVGF6KGXJ/m1Up1AdenEY6OiBmVVrR00GgmkjqrzUYoItenebVSXcD1aTbunjIzs9wcNMzMLDcHjbFzbaMLUGOuT/NqpbqA69NUPKZhZma5uaVhZma5OWiYmVluDhojJOkoSQ9IekbSSkl/mNIPk/RdSavS30NTuiR9TVK3pKcknViS19lp+1WSzq52zDrX50BJj0p6MtXnL1P6HEmPpHLflm5lT7rd/W0p/RFJs0vyuiSld0k6oxH1KSnLOEkrJH0zPS9sfSQ9L+lpSU9I6kxpRX2/TZV0h6SfSHpW0ikFrstx6X8y8HhF0h8VtT5Digg/RvAAZgInpuVDgOeAecD/AC5O6RcDl6fls4BvAQJOBh5J6YcBq9PfQ9PyoQ2oj4DJaXkC8Egq5+3AopR+DfCZtPxZ4Jq0vAi4LS3PA54EJgFzgJ8C4xr4f7oIuAX4Znpe2PoAzwPTy9KK+n67ETgvLU8Epha1LmX1Ggf0AEe3Qn0q1rHRBWiVB9mPRX0Q6AJmprSZQFda/gdgccn2XWn9YuAfStL3265BdXkT8CPg3WRXro5P6acA96Xl+4BT0vL4tJ2AS4BLSvLat10D6jEL+B7wAeCbqXxFrs/zvDFoFO79Rvarm2tIE3GKXJcKdTsd+H6r1KfSw91TNZC6Mt5FdnZ+eERsSKt6gMPT8pHAupLd1qe0auljLnXlPAG8BHyX7Kx6S0T0VSjbvnKn9VuBaTRRfYCvAn8G9Kfn0yh2fQL4jqTHJZ2f0or4fpsD9AL/J3UdXifpYIpZl3KLgFvTcivU5w0cNEZJ0mTgX4E/iohXStdFdrpQmDnNEbE3Ik4gO0M/CXhbg4s0YpJ+C3gpIh5vdFlq6D0RcSJwJnCBpPeWrizQ+208cCJwdUS8C9hO1n2zT4Hqsk8aH/sQ8C/l64pYn2ocNEZB0gSygHFzRCxNyS9KmpnWzyQ7awd4ATiqZPdZKa1aesNExBay32o/BZgqaeBngUvLtq/caf0UYCPNU59fAz4k6XlgCVkX1f+muPUhIl5If18C7iQL7EV8v60H1kfEI+n5HWRBpIh1KXUm8KOIeDE9L3p9KnLQGCFJIvvd8mcj4u9KVi0DBmY9nE021jGQ/ok0c+JkYGtqut4HnC7p0DS74vSUNqYkzZA0NS0fRDY+8yxZ8FiYNiuvz0A9FwL3p7OpZcCiNBtpDjAXeHRsavG6iLgkImZFxGyyLoP7I+JjFLQ+kg6WdMjAMtn75McU8P0WET3AOknHpaRTgWcoYF3KLOb1rikofn0qa/SgSlEfwHvImptPAU+kx1lk/eDfA1YB/w84LG0v4EqycYKngY6SvD4FdKfHJxtUn3cAK1J9fgx8MaUfQ/Yl2U3W7J6U0g9Mz7vT+mNK8vpCqmcXcGYT/K/ez+uzpwpZn1TuJ9NjJfCFlF7U99sJQGd6v91FNluokHVJ5TiYrGU6pSStsPUZ7OHbiJiZWW7unjIzs9wcNMzMLDcHDTMzy81Bw8zMcnPQMDOz3Bw0zMwsNwcNMzPL7f8DDv5Utjzsj74AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Baseline\n",
    "wuhan_case_Jan21 = Estimate_wuhan_case(model_name = 'Baseline', date = '2019-01-21', \n",
    "                                      wuhan = Wuhan(population=11000000, \n",
    "                                                    airportCatchment=19000000, \n",
    "                                                    internationalTraveller=3301),\n",
    "                                      international = International(cases = 7),\n",
    "                                      coronavirus = Coronavirus(incubation=6, \n",
    "                                                                onsetTodetection=4))\n",
    "\n",
    "print(wuhan_case_Jan21)\n",
    "# Plot the distrubution of estimated Coronavirus cases in Wuhan \n",
    "wuhan_case_Jan21.plot_distribution()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sensitivity Analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2019-01-21 \n",
      "Scenario: Baseline\n",
      "Exported number of confirmed cases: 7\n",
      "Detection Window: 10\n",
      "Time from onset of symptoms to detection: 4\n",
      "Daily International passengers travelling out of Wuhan: 3301\n",
      "Effective catchment population of Wuhan Interntional Airport: 19000000\n",
      "\n",
      "Estimated number of cases in Wuhan: 4029.0820963344436\n",
      "Estimated 95% confidence interval: (1730, 7780)\n",
      "\n",
      "------------------------------------------------------------------------------------\n",
      "\n",
      "2019-01-21 \n",
      "Scenario: Smaller catchment\n",
      "Exported number of confirmed cases: 7\n",
      "Detection Window: 10\n",
      "Time from onset of symptoms to detection: 4\n",
      "Daily International passengers travelling out of Wuhan: 3301\n",
      "Effective catchment population of Wuhan Interntional Airport: 11000000\n",
      "\n",
      "Estimated number of cases in Wuhan: 2332.6264768252045\n",
      "Estimated 95% confidence interval: (1000, 4490)\n",
      "\n",
      "------------------------------------------------------------------------------------\n",
      "\n",
      "2019-01-21 \n",
      "Scenario: Shorter detection window\n",
      "Exported number of confirmed cases: 7\n",
      "Detection Window: 8\n",
      "Time from onset of symptoms to detection: 2\n",
      "Daily International passengers travelling out of Wuhan: 3301\n",
      "Effective catchment population of Wuhan Interntional Airport: 19000000\n",
      "\n",
      "Estimated number of cases in Wuhan: 5036.352620418054\n",
      "Estimated 95% confidence interval: (2160, 9720)\n",
      "\n",
      "------------------------------------------------------------------------------------\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# 2019-01-21: same result as papers \n",
    "sensitivity_analysis(date='2019-01-21', wuhan_population=11000000,\n",
    "                     airportCatchment=19000000,international_case=7, \n",
    "                     onsetTodetection=4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2019-01-22 \n",
      "Scenario: Baseline\n",
      "Exported number of confirmed cases: 8\n",
      "Detection Window: 10\n",
      "Time from onset of symptoms to detection: 4\n",
      "Daily International passengers travelling out of Wuhan: 3301\n",
      "Effective catchment population of Wuhan Interntional Airport: 19000000\n",
      "\n",
      "Estimated number of cases in Wuhan: 4604.6652529536495\n",
      "Estimated 95% confidence interval: (2100, 8550)\n",
      "\n",
      "------------------------------------------------------------------------------------\n",
      "\n",
      "2019-01-22 \n",
      "Scenario: Smaller catchment\n",
      "Exported number of confirmed cases: 8\n",
      "Detection Window: 10\n",
      "Time from onset of symptoms to detection: 4\n",
      "Daily International passengers travelling out of Wuhan: 3301\n",
      "Effective catchment population of Wuhan Interntional Airport: 11000000\n",
      "\n",
      "Estimated number of cases in Wuhan: 2665.8588306573765\n",
      "Estimated 95% confidence interval: (1220, 4940)\n",
      "\n",
      "------------------------------------------------------------------------------------\n",
      "\n",
      "2019-01-22 \n",
      "Scenario: Shorter detection window\n",
      "Exported number of confirmed cases: 8\n",
      "Detection Window: 8\n",
      "Time from onset of symptoms to detection: 2\n",
      "Daily International passengers travelling out of Wuhan: 3301\n",
      "Effective catchment population of Wuhan Interntional Airport: 19000000\n",
      "\n",
      "Estimated number of cases in Wuhan: 5755.831566192062\n",
      "Estimated 95% confidence interval: (2630, 10700)\n",
      "\n",
      "------------------------------------------------------------------------------------\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# 8 exported cases until Jan 22 \n",
    "# 29 exported cases until Jan 26 \n",
    "sensitivity_analysis(date='2019-01-22', wuhan_population=11000000,\n",
    "                        airportCatchment=19000000,international_case=8, onsetTodetection=4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2019-01-29 \n",
      "Scenario: Baseline\n",
      "Exported number of confirmed cases: 67\n",
      "Detection Window: 10\n",
      "Time from onset of symptoms to detection: 4\n",
      "Daily International passengers travelling out of Wuhan: 3301\n",
      "Effective catchment population of Wuhan Interntional Airport: 19000000\n",
      "\n",
      "Estimated number of cases in Wuhan: 38564.071493486816\n",
      "Estimated 95% confidence interval: (30000, 48470)\n",
      "\n",
      "------------------------------------------------------------------------------------\n",
      "\n",
      "2019-01-29 \n",
      "Scenario: Smaller catchment\n",
      "Exported number of confirmed cases: 67\n",
      "Detection Window: 10\n",
      "Time from onset of symptoms to detection: 4\n",
      "Daily International passengers travelling out of Wuhan: 3301\n",
      "Effective catchment population of Wuhan Interntional Airport: 11000000\n",
      "\n",
      "Estimated number of cases in Wuhan: 22326.56770675553\n",
      "Estimated 95% confidence interval: (17340, 28030)\n",
      "\n",
      "------------------------------------------------------------------------------------\n",
      "\n",
      "2019-01-29 \n",
      "Scenario: Shorter detection window\n",
      "Exported number of confirmed cases: 67\n",
      "Detection Window: 8\n",
      "Time from onset of symptoms to detection: 2\n",
      "Daily International passengers travelling out of Wuhan: 3301\n",
      "Effective catchment population of Wuhan Interntional Airport: 19000000\n",
      "\n",
      "Estimated number of cases in Wuhan: 48205.08936685853\n",
      "Estimated 95% confidence interval: (37510, 60600)\n",
      "\n",
      "------------------------------------------------------------------------------------\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# 2019-01-29: 67 confirmed cases overseas\n",
    "sensitivity_analysis(date='2019-01-29', wuhan_population=11000000,\n",
    "                     airportCatchment=19000000,international_case=67, \n",
    "                     onsetTodetection=4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
